function [An,Dn] = InclinationChange(Coe,u,di,Units)
    
    a = Coe(1);
    %[km]Semimajor axis.
    
    e = Coe(2);
    %[]Eccentricity.
    
    Om = Coe(4) * pi / 180;
    %[rad]Right ascension of the ascending node.
    
    w = Coe(5) * pi / 180;
    %[rad]Argument of periapsis.
    
    p = a * (1 - e^2);
    %[km]Semiparameter.
    
    n = sqrt(u / a^3);
    %[rad/s]Mean angular speed.
    
    theta(1) = Coe(6) * pi / 180;
    %[rad]Initial satellite true anomaly.
    
    theta(2) = -w;
    %[rad]True anomaly of the AN.
    
    theta(3) = pi - w;
    %[rad]True anomaly of the DN.
    
    E(1) = 2 * atan(sqrt((1 - e) / (1 + e)) * tan(theta(1) / 2));
    %[rad]Initial satellite eccentric anomaly.
    
    E(2) = 2 * atan(sqrt((1 - e) / (1 + e)) * tan(theta(2) / 2));
    %[rad]Eccentric anomaly of the AN.
    
    E(3) = 2 * atan(sqrt((1 - e) / (1 + e)) * tan(theta(3) / 2));
    %[rad]Eccentric anomaly of the DN.
    
    t(1) = (E(1) - e * sin(E(1))) / n;
    %[s]Initial satellite time since periapsis passage.
    
    t(2) = (E(2) - e * sin(E(2))) / n;
    %[s]Time since periapsis passage of the AN.
    
    t(3) = (E(3) - e * sin(E(3))) / n;
    %[s]Time since periapsis passage of the DN.
    
    N = [cos(Om); sin(Om); 0];
    %[]Line of nodes direction in ECI coordinates.
    
    R = Rotation(N,di,Units);
    %[]Matrix that rotates the velocity vector during the inclination change.
    
    Tpqw2eci = Pqw2Eci(Coe);
    %[]Matrix that transforms vectors from PQW coordinates to ECI coordinates.
    
    An.Vminus = sqrt(u / p) * Tpqw2eci * [-sin(theta(2)); cos(theta(2)) + e; 0];
    %[km/s]AN velocity WRT the Earth in ECI coordinates before the inclination change.
    
    An.Vplus = R * An.Vminus;
    %[km/s]AN velocity WRT the Earth in ECI coordinates after the inclination change.
    
    An.dV = An.Vplus - An.Vminus;
    %[km/s]AN inclination change delta-v vector in ECI coordinates.
    
    An.dv = norm(An.dV);
    %[km/s]AN inlcination change delta-v magnitude.
    
    An.dt = t(2) - t(1);
    %[s]Time that it takes the satellite to reach the AN.
    
    if An.dt < 0
        
        T = 2 * pi / n;
        %[s]Orbital period.
        
        An.dt = An.dt + T;
        %[s]Makes the time of flight positive if needed.
        
    end
    
    Dn.Vminus = sqrt(u / p) * Tpqw2eci * [-sin(theta(3)); cos(theta(3)) + e; 0];
    %[km/s]DN velocity WRT the Earth in ECI coordinates before the inclination change.
    
    Dn.Vplus = R * Dn.Vminus;
    %[km/s]DN velocity WRT the Earth in ECI coordinates after the inclination change.
    
    Dn.dV = Dn.Vplus - Dn.Vminus;
    %[km/s]DN inclination change delta-v vector in ECI coordinates.
    
    Dn.dv = norm(Dn.dV);
    %[km/s]DN inlcination change delta-v magnitude.
    
    Dn.dt = t(3) - t(1);
    %[s]Time that it takes the satellite to reach the DN.
    
    if Dn.dt < 0
        
        T = 2 * pi / n;
        %[s]Orbital period.
        
        Dn.dt = Dn.dt + T;
        %[s]Makes the time of flight positive if needed.
        
    end
    
end
%===================================================================================================